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ALTERNATIVE TO EVOLVING SURFACE FINITE ELEMENT 

METHOD 

MARYIA BORUKHAVA * AND HEIKO KRONERt 

Abstract. ESFEM is a method introduced in [5] in order to solve a linear advection-diffusion 
equation on an evolving two-dimensional surface with finite elements by using a moving grid with 
nodes sitting on and evolving with the surface. The evolution of the surface is assumed to be given 
as a smooth one-parameter family of embeddings of a fixed initial surface into R 3 satisfying uniform 
C 4 bounds. We calculate an equivalent transformed equation which is defined on the fixed initial 
surface and can hence be solved numerically on a fixed grid. We present numerical examples which 
indicate that both approaches are essentially of the same accuracy. 
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1. Introduction. In many applications it is important to consider PDEs which 
are defined on surfaces and not in Euclidean space, especially in the case of parabolic 
equations it is of interest to assume that these surfaces (where the equation is defined) 
evolve with respect to time in a certain prescribed way. In [5] the so-called evolving 
surface finite element method (ESFEM) is proposed in order to solve an advection- 
diffusion equation on an evolving surface, cf. i5j Sections 1.1 and 1.2]. This setting 
models e.g. the transport of an insoluble surfactant on the interface between two 
flowing fluids or pattern formation on the surfaces of growing organisms modeled 
by reaction-diffusion equations, cf. 0 Section 1.4] for further and a more detailed 
exposition of applications. There are several papers which deal with linear parabolic 
equations on evolving surfaces, e.g. in mm it is shown that classical L 2 - and L°°- 
estimates carry over to ESFEM and in mm it is shown that this also holds for error 
estimates of Runge-Kutta schemes and Backward difference schemes; we also mention 

mm- 

We present the idea of ESFEM according to 0 and note that this method is intro¬ 
duced therein in order to solve a special linear parabolic equation, namely Equation 
OD, with the finite element method. 

Let r 0 be a smooth, closed, connected and oriented hypersurface in R 3 . Let 
<h(f, •) : r 0 -A R 3 , t £ [0,T 0 ], T 0 > 0, be a family of embeddings, $ smooth, T{t) = 
<h(f, r 0 ) the moving surfaces. We define the set 

(1.1) g To = 1J Wxr(t) 

te[o,T 0 ] 

and consider there the advection-diffusion equation 

(1.2) u + uV r(t) • v - V r(t) • (T> 0 V r(t) u) = 0 

where u denotes the material derivative, v(p t ,t) = ^|(<h(f,-) _1 (p t ), t) the velocity of 
the moving surface in a point ( t,pt ) £ Gt 0 and hence V r ^ -v its tangential divergence. 
Furthermore, the diffusion coefficient Dq : Gt 0 —> Mat(n+ 1,R) is so that it vanishes 
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on the normal space of the moving surfaces, i.e. -D 0 (t,p t );/ = 0 f° r a ll v G N pt where 
N pt is the normal space of T(i) in p t . Equation (11.21) can be written in variational 
form as 

(1.3) 4/ up + f A,V r(t) • V r (‘V = f up Vp e C°°(G To )- 

Jr(t) Jr(t) Jr(t) 

ESFEM considers at every time t an interpolating polyhedral surface which 

approximates the evolving surface T(t) and which consists of triangles with vertices 
Xj(t) = $(t, Xj(0)), j = 1 sitting on T(f) and moving with the surface. Here, 

Xj(Q), j = 1 are fixed nodes on the initial surface T(0). The finite element 

basis functions pi , i = are defined on 

(1-4) G h To = U (t}xr k (f) 

tG [0,Tb] 


and chosen so that pi(t , •) is piecewise linear (i.e. linear on each triangle of Th(t)) 
with pi(t,Xj(t)) = Sij. 

The fully discrete scheme from [5] in order to solve (O which uses ESFEM can 
be found in [5] Equation (7.2)] and is as follows. Let to = 0 < ... < £m = To, M £ N, 
be a partition of the time interval, T™ = Th{t m ), u m = u(-,t m ) for a function u on 
Gy o and u ° € T/j(0) an initial function (and for simplicity Dq the identity) then we 
solve for m = 0,..., M — 1 the linear system 


(1.5) 


< + vr +1 


+ [ vC +1 u ™+i. v c + V7 +1 

J T m + l 


=1/ 
r J rr 


<P?, j = h-,N. 


Our different approach to solve (11.21) with the finite element method reformulates 
Equation (11.21) as an equivalent equation on the surface T(0), cf. Equation (13.131) . 
and uses a fixed grid (the one consisting of the nodes X 7 (0), j = 1, ...,1V,) for the 
finite element approximation (again an implicit Euler method) which leads to a finite 
element solution ft™ defined on T^(0). Obviously, both approaches coincide if $(t, •) 
is the identity, i.e. if there is no motion. 

Our aim is to use both approaches, i.e ESFEM and the method which uses the 
transformed equation, in some example cases and to provide the corresponding error 
tables. 

The paper is organized as follows. In Section [2] we recall some facts about hy¬ 
persurfaces in the Euclidean space, in Section [3] we derive the reformulated equation, 
in Section [4] we present our chosen examples and a more detailed description of the 
parts from which the coefficients of the transformed equation are put together in the 
implementation. The error tables of the numerical calculations can be found in Tables 
1 to 6. 


2. Hypersurfaces in K"+ 1 . We recall some facts and notations of embedded 
hypersurfaces in R ra+1 from |7]. Let F : Q —> R n+1 with C K n open be a smooth 
embedding and M = F(Cl). For p £ H the coordinate tangent vectors diF(p) = (p), 

1 < i < n, provide a basis of the tangent space T X M at x = F(p). 

The metric on M is given by 


(2.1) 


9ij = diF ■ djF 
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for 1 < i,j < n, the inverse metric by ( g 13 ) = (g,y) 1 . 

The tangential gradient of a function h : M —> M is defined by 

( 2 . 2 ) X M h = g i3 djhdiF 
where we sum over repeated indices. 

For a smooth tangent vector field X = X l diF = g 13 XjdiF on M (note that 
Xi = X ■ diF) we define the covariant derivative tensor by 

( 2 . 3 ) VfX 3 = d z X 3 + T\ k X k 
where the Christoffel symbols r* • are given by 

(2-4) r£- = ^g kl (digji + djga - digtj). 

The tangential divergence of X on M is defined by 
( 2 . 5 ) V M • X = div M X = XfX i 


and the Laplace-Beltrami operator of h on M by 

(2.6) A M h = div M V M /i = V M • (V M /i). 

For a smooth vector field X : M —> R" +1 which is not necessarily tangent on M 
we can also define the divergence with respect to M by 

(2.7) V M • A = div M X = g i3 diX ■ d 0 F 


which reduces to the above expression if X is tangent on M. We remark that, of 
course, the divergence in (12.71) is (and transforms like) a scalar function (when chang¬ 
ing local coordinates of M). Furthermore, the tangential gradient transforms like a 
scalar function when changing coordinates in M. 

3. Reformulation of the equation on a fixed surface. In this section we 
derive Equation (13.131) which is an equivalent reformulation of (11.21) on T(0). The 
calculation is straightforward but we still present details. 

Let fl, VL C R” be open and $ : a diffeomorphism. The linear differential 

operator L : iJ 2 (fi) —> L 2 (fl) 

(3.1) Lu = a 13 DiDjU + b l DiU + cu 

with coefficients a 13 ,b l ,c € L°°(f2), a 13 symmetric, transforms into the operator L : 
H 2 (n) -S- L 2 (n) 

(3.2) Lu = a 13 DiDjU + b'Dj/i + cu, 

i.e. we have (Lu) o $ = Lu for the quantity u = u o <h. Here, by writing x(x) = d>(x) 
and x(x) = $” 1 (x) we set 


(3.3) 


V3 =n rs 


b k =b n 


dx l dx 3 
dx r dx s 

dx k ,, d 2 x m 
- a 13 - 


dx n 


dx r dx s dx k 


c =c 


dx r dx s dx 1 dxi dx 171 






4 


A remark on the evolving surface finite element method 


and o$ operations are suppressed. 

These formulas carry over to the surface case when using local coordinates. Let 
now fl, fl be (open subsets of) hypersurfaces in R 3 , $ : Cl —> fl a diffeomorphism and 
let a 1 - 3 , b 1 and c be L°°-sections of the tensor bundles T 2 ’°(fl), T 1,0 (fl) and T 0,0 (fl), 
respectively, and assume that a v is symmetric. Let L be defined by 


Lu = a IJ Y''Y''u + b‘V?u + cu, u G iL 2 (fl), 


(3.4) 


where V s2 denotes the Levi-Cevita connection with Christoffel symbols T*- on fl (and 
V n and T% correspondingly on Cl) then we have 


Lu =a l i DiDjU + ( b k — a l ^T k -)DkU + cu 
=a*- J DiDjU + b k D k u + cu 


(3.5) 


where Di denote ordinary partial derivatives. Using the formulas (13.31) we get a 
transformed operator 


LU = a ij C7?C7fu + VV?u + cu, UG H 2 {Cl) 


(3.6) 


where now in local coordinates (x l ) of fl and (x l ) of fl we have 



rk Trn dx k d 2 x m dx r dx s dx k 

b k = Jfri -- 


(3.7) 


dx m dx r dx s dx l dxi dx m 



c =c. 


A choice of local coordinates in Cl induces via 4> local coordinates in fl and we stipulate 
that in the following the local coordinates of fl and fl are related in this way. Then 
the formulas (13.71) simplify to 


(3.8) 


b k =b k + a^f k J 


c =c. 


Let us consider the case where the main part is in divergence form 


(3.9) 


Lu =Vf(a ij Vpu) + b‘Vfu + cu 

+ (\/fa ij + 6*)vf u + cu 


then we get the transformed operator with main part in divergence form 


Lu —+ b‘V?u + cu 

=Vp {a ij V? u) + (b* - Vpa ij )vf u + cu 
=Vf(a i3 vju) + b‘wfu + cu 


(3.10) 
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where 

a ij =a ij 

(3.11) b k =b k + a i3 T*) + V)V fc = b k + a y (f£ - T k 0 ) + V?a ifc 

C —C 

and hence 

(3.12) b k =b k + a iJ { f£- - I*) + V)V fc - Vfa ik . 

Let us assume that u and $ are as in Section [T] and that D 0 is the identity. We 
define u(t,x) = it(f, 4>(t, x)) and transform (11.21) in the equivalent equation 

^-u- v[ (0 V J (f)Vy (0) fi) 

(3.13) dt 5 

+ (ff y (i)(T£.(i) - r^-(O)) + V?V(i))V?°>fi + fiv r « • « = 0 

on [0, To] x r(0) where the time derivative is now a usual partial derivative, gij{t) 
denotes the metric and iW (t) the Christoffel symbols of r(f). Here, we used that the 
covariant derivative of the metric vanishes and the coupling of local coordinates via 

Remark 3.1. We note that if one considers instead of (El a general linear 
parabolic equation 

(3.14) u - a ij v[ (0 Vj (t) w + 1[t) u + cu = f 

on Gt 0 where 

(3.15) : G To ->• T 2 '°(r(f)), V : G To T 1 - 0 ^)), c, / : G To -»• K 

so that a 1 - 7 (f, •) and b l (t , •) are sections of the corresponding bundles then the transfor¬ 
mation rules (13.111) and (13.121) - of course - provide a reformulation of this equation 
on the fixed surface as well. 

4. Examples and implementation. We choose T(0) = d£i(0) C R 3 , T 0 = § 
and define the motion of the surface by 


(4.1) $(t,x)=A(t)x 
where we consider the cases 

/ \/l + 5sinf 0 0\ 

(4.2) A(t) = 0 10, 

V 0 0 1/ 

/l + txf 0 0\ 

(4.3) A(t) = 0 10 

V 0 0 1/ 

and 

( cos(r}n(x 3 )t) - sin(?^(x 3 )f) 0\ 

sm(r/g(x 3 )t) cos(rgi(x 3 )t) 0 

0 0 l) 
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with 77 = 5 and ^(£ 3 ) = ^£ 3 . We append a right-hand side to the equation (11.31) 
(and correspondingly to (13 131) ') and set an initial value which are so that the exact 
solution of the transformed equation is in all cases u(t,x) = e~ 6t X\X 2 ■ We calculated 
the parts from which the three right-hand sides (corresponding to the three cases) 
can be put together by hand and put them together within the implementation. For 
the discretization we use the coupling At = A£ 2 between the step sizes in time and 
space. 

We point out that the third example is chosen to observe the phenomenon of 
having only tangential motion which deteriorates the mesh in ESFEM (in the sense 
that the ratio of the diameter and the incircle radius of a triangle might become large) 
and of course also affects the coefficients in the transformed equation. We remark that 
for ESFEM the (relative) mesh size may change and tangential motion of the nodes 
may deteriorate the mesh. While the former might be compensated by appending 
additional nodes the latter can be compensated by introducing additional tangential 
motion of the mesh (ALE-ESFEM), see [9] for the latter. 

To be able to input the transformed coefficients as stated in (13.131) into the Dis¬ 
tributed and Unified Numerics Environment (DUNE), see [TJ[2J[3], we present some 
auxiliary facts. Let x = (x l )i<i <3 be Euclidean coordinates in R 3 and £ = (£ J )i<j <2 
local coordinates in T(0). Then x = £*(£-') can be seen as a local representation of 
the embedding of T(0) into R 3 and <F(£*(£ J ), t) of the embedding of T(t) into R 3 . For 
the latter we drop the time dependence for simplicity and write $ o x = $(£ l (^)). 
We have 


(4.5) 


d d 




AL_(#o*).A(4.o I ) + 


d 2 


d 

(<f> O x) • —-( 4 > O X), 


d^d^ k v ' d£ 


furthermore, 

(4.6) V£(«V(i) = -g ik (t)g lj (t)d r g kl (t) + T(0)Ug li + r(0)j jff K . 


We calculate the tangential divergence of the evolution speed v 


A(t)x 


( 4 -7) V r(t) • v = g lJ (t)(-^A(t)x + A(t)-^x) ■ -^-($ o x). 

Let £ = ( x l ) C T(O) be given. Let a : {1,2,3} — > {1,2,3} be a bijection so that 
|£ CK ( 1 ) | = mirij |x“W|, £°4 3 ) ^ 0 and set x = x a ^ 3 \ y = x a ( 2 \ z = x a ^ , ip = £ 2 = 
arctan(^) and 


(4.8) 


o=e 


arccos (z) if x > 0 
27t — arccos(^) else. 


Then we have 


(4.9) 


£ = sin 0 cos ip 
y = sin 6 sin ip 
z = cos 6 
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and 


(4.10) 


dx 

w 

dx 

w 

dy 

dC 

dy 

dC 


= cos 9 cos p 
= — sin 9 sin ip 
= cos 9 sin ip 
= sin 9 cos ip 


dz 

w 

dz 

W 


= — sin 9 


= 0 . 


We have to stipulate what we understand for a given triangle T of r^O) by a repre¬ 
sentation of a*- 7 and b l in T (where T denotes the lift of T to T(0)) with respect to 
given orthonormal basis v\ , V 2 in T. For this purpose we let v\ , V 2 be the orthogonal 
projections of jjprX and -gjrzX onto the plane which contains T. Then we define /?(., 
k, l = 1 , 2 by 


vi =0[v i + 0 lv 2 
V2 =02^ 1 + I3%V2 


then 


(4.12) 


~b l = b^3), a ij = a kl Pi/3j 


are what we understand by evaluating a lJ ,b l with respect to vi,V 2 - Let P denote 
the plane containing T and let M be the matrix representation with respect to the 
standard basis in R 3 of an (arbitrary) extension to M 3 x R 3 of the bilinear form on 
P x P represented by the matrix a y with respect to the basis v±,V 2 - Correspondingly 
we construct a vector B £ M 3 from b l . Then M and B are the coefficients which are 
compatible with the input format of DUNE. 

In order to calculate the right-hand sides of the transformed equation for the 
exact solution u (this means we evaluate the left-hand side of (13.131) for our special 
u) we use our coordinates (£*) to calculate the summand 

(4.13) ff ij (t)v[ (0) Vj (0) M 

which can be put together from 


(4.14) 


vr ( °\ r(o) «= 


d 2 


d^d^ 




Note, that u = e~ 6t XiX 2 ^ e~ 6t xy in general (depending on a). 

We use the mesh generator Gmsh, see [TO] , and for the implementation, more 
precisely, DUNE-FEM, a discretization module for solving PDEs which depends on 
DUNE. 
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